---
title: "R Script: Authoritarianism, anti-egalitarianism or economic conservatism? The role of social and economic attitudes in supporting environmental protection in the United Kingdom"
name: Pablo Sallabera & Araceli Mateos
output: html_document
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


# Basic information
 Author:    Pablo Sallabera & Araceli Mateos
 Contact:   sallabera@gmail.com
 File name: sallabera_mateos_2022_RSript.rmd
 Title:   "Authoritarianism, anti-egalitarianism or economic conservatism? The role of social and economic attitudes in supporting environmental protection in the United Kingdom"
 Started:   2021-01-02
 Summary:   RMD R-Script

# Packages
```{r}

# To install the packages use these functions:

# install.packages("tidyverse")
# install.packages("haven")
# install.packages("dplyr")
# install.packages("forcats")
# install.packages("ggplot2")
# install.packages("skimr")
# install.packages("sjlabelled")
# install.packages("readxl")
# install.packages("huxtable")
# install.packages("stats")
# install.packages("gmodels")
# install.packages("tableone")
# install.packages("kableExtra")
# install.packages("margins")
# install.packages("Hmisc")
# install.packages("janitor")
# install.packages("purrr")
# install.packages("naniar")
# install.packages("summarytools")

# install.packages("sjPlot")
# install.packages("texreg")
# install.packages("foreign")
# install.packages("jtools")
# install.packages("texreg")
# install.packages("rms")
# install.packages("stats")
# install.packages("magrittr")
# install.packages("nnet")
# install.packages("stargazer")
# install.packages("oglmx")
# install.packages("fmsb")
# install.packages("ggpubr")
# install.packages("Hmisc")
# install.packages("psych")
# install.packages("broom.mixed")




# For loading the packages use these others functions: 

library(tidyverse)
library(haven)
library(forcats)
library(dplyr)
library(ggplot2)
library(skimr)
library(sjlabelled)
library(readxl)
library(huxtable)
library(stats)
library(gmodels)
library(tableone)
library(kableExtra)
library(margins)
require(Hmisc)
library(janitor)
library(purrr)
library(naniar)
library(summarytools)

library(sjPlot)
library(texreg)
library(foreign)
library(jtools)
library(texreg)
library(rms)
library(stats)
library(magrittr)
library(nnet)
library(stargazer)
library(oglmx)
library(fmsb)
library(ggpubr)
library(psych)
library(broom.mixed)


```


# Source full data 
```{r}

# Here is the full dataset with all the varaibles and waves. In the next step we reduce it to make it more manageable.you can find the reduced dataset in the working directory ("data_bes_red.sav").


  # bes_full <-read_sav("BES2019_W20_Panel_v0.1")
# We mutate this function because takes a long time to run. We recomend to use the "data_bes_red.sav" dataset that is shown below.  
```



# Selection of varaibles. 

Note: you can download the questionnaire here in order to have more information about the variables and their labels. For the majority of the variables we use the same names: https://www.britishelectionstudy.com/wp-content/uploads/2020/11/Bes_wave20Documentation_V2.pdf

```{r}

# We select the varaibles of interest (this function is muted):

 # df_bes_full <- bes_full %>% select(id, 
                   wave7, # variables that indicate if the interviewee has answered that wave of the survey. 
                   wave14,
                   wave17,                          
                   wt_new_W7, # weighting coefficient
 
                ## Dependent variables 
                   enviroGrowthW7, # Economy VS Environment scale W1
                   enviroGrowthW17, # Economy VS Environment scale W2
                   enviroProtectionW7, # Environment Protection W1
                   enviroProtectionW14, # Environment Protection W2
                   
                ## Sociodemographics
                   ageW7,  
                   gender,
                   p_educationW7,
                   
                ## Socioeconomic Status
                   p_work_statW7, # Employment Status 
                   econPersonalRetroW7, # assessment of the evolution of the personal economic situation
                   econGenRetroW7, # subjective assessment of the evolution of the economic situation of the country 
                   p_gross_householdW7, # household income
                   ns_sec_analyticW6W7W8W9, # occupation 
                   subjClassW2_W4W7W9, # subjective social class               
                
                ## Media consumption 
                   infoSourceTVW7,
                   infoSourcePaperW7,
                   infoSourceRadioW7,
                   infoSourceInternetW7,
                
                ## Party Identity
                   partyIdW7, 
                   leftRightW7, # ideology  

 ## Heath- Evans: Authoritarianis items 
    lr1W7_W9, # Government should redistribute 
    lr2W7_W9, # Big business takes advantage of ordinary people
    lr3W7_W9, # Ordinary working people do not get their fair share of the nation’s wealth
    lr4W7_W9, #There is one law for the rich and one for the poor
    lr5W7_W9, # Management will always try to get the better of employees if it gets the chance
 ## Heath- Evans: Economic Conservatism items 
    al1W7_W9, # Young people today don’t have enough respect for traditional British Attitudes
    al2W7_W9, # For some crimes, the death penalty is the most appropriate sentence
    al3W7_W9, # Schools should teach children to obey authority 
    al4W7_W9, # Censorship of films and magazines is necessary to uphold moral standards
    al5W7_W9, # People who break the law should be given stiffer sentences
            
                
   ## Anti-Egalitarian items
      blackEqualityW6_W14, # Attempts to give equal opportunities to ethnic minorities
      femaleEqualityW6_W14, # Attempts to give equal opportunities to women
      gayEqualityW6_W14, # Attempts to give equal opportunities to gays and lesbians
                
   ## Services (2nd operationalisation): Economic Conservatism items
      cutsTooFarNationalW7, #Cuts to public spending in general
      cutsTooFarNHSW7, #Cuts to NHS spending*
      cutsTooFarLocalW7, #Cuts to local services in my area
      privatTooFarW7, #Private companies running public services
            
                ) %>% filter(wave7==1) # We delete those interviewees that have not answered the wave7 to reduce the dataset size.   
 


# We start asignating the label NA to those Attitudes that we want to eliminate. In the Party Identity variable we eliminate the minor and regional parties:  
        df_bes_full$partyIdW7[df_bes_full$partyIdW7==5] <- NA
        df_bes_full$partyIdW7[df_bes_full$partyIdW7==8] <- NA
        df_bes_full$partyIdW7[df_bes_full$partyIdW7==11] <- NA
        df_bes_full$partyIdW7[df_bes_full$partyIdW7==9] <- NA

      
# We run this function that assign NA easily, but this function takes several minutes: 
     # df_bes_red <- df_bes_full %>% replace_with_na_all(condition = ~.x %in% c(9999, 998, 999))%>%
    replace_with_na(replace = list(p_educationW7    = c(19, 20),
                                   p_gross_householdW7    = c(16, 17))) 


```

# Reduced Data Set 
```{r}

# This is the reduced data set that you can find in the working directory: 
  # write_sav(df_bes_red, "data_bes_red.sav")  

# Here to upload the reduced data set: 
df_bes <-read_sav("data_bes_red.sav")


# With this function you can preview the descriptive statistics of the data set. 
 view(dfSummary(df_bes))
```



# 1. LABELING  

Note: You can use the drop-down menu on the right to see all the sections, organized in an index. 

## 1.1. Independent Variables 
```{r}

# Function for assign labels to the factor varaibles: 
 df_bes <- df_bes %>%  dplyr::mutate(
  gender=as_label (gender),
  p_work_statW7=as_label (p_work_statW7),
  p_work_statW7=fct_collapse(p_work_statW7, 
    "Unemployed"=c("Not working", "Unemployed"), # We join several redundant Attitudes
    "Working"=c("Working full time (30 or more hours per week)", "Working part time (Less than 8 hours a week)", "Working part time (8-29 hours a week)" )), 
                
  partyIdW7=as_label(partyIdW7)) 
                                        
                    
```


## 1.2. labeling Subcategories
```{r}
# Occupation

# These are the cathegories you can find in the questionnaire: 

# 11 ○ Employers in large organisations and higher managerial
# 12 ○ Higher professional occupations
# 20 ○ Lower professional and managerial and higher supervisory
# 30 ○ Intermediate occupations
# 40 ○ Employers in small organisations and own account workers
# 50 ○ Lower supervisory and technical occupations
# 60 ○ Semi-routine occupations
# 70 ○ Routine occupations

# We join some of these groups in 5 main cathegories: 
df_bes <-df_bes %>% 
       mutate(occupation=fct_collapse(as.factor(ns_sec_analyticW6W7W8W9), 
                             "High-quali"=c("11", "12"), 
                             "Average-quali 1"=c(20),
                             "Average-quali 2"=c(30),                                                                                           "freelancers/small-business"=c(40),
                             "Low-quali"=c(50, 60, 70)))

# Use this function to see the distribibution of each cathegory:
table(df_bes$occupation)





# We rellabel some Attitudes of the Subjective Social Class varaible  
df_bes <-df_bes %>% 
       mutate(subjClassW7.f=fct_collapse(as.factor(subjClassW2_W4W7W9), 
            "Middle-class"=c("1"), 
            "Working-class"=c("2"), 
            "Other-class"= c("3"),
            "No-class"= c("0")))

table(df_bes$subjClassW7.f)




# We proceed in the same way with the education variable:
 df_bes <-df_bes %>% mutate(
           educationW7= fct_collapse(as.factor(p_educationW7), 
           "No university"= c(1:14),
           "University/High qualific"= c("15", "16", "17"), 
           "Other high quali" =c(18)))
 
table(df_bes$educationW7)
  



     
# Party Identity relabeling: 

# We change the label name of "No - none" to "No-partyId". 
   df_bes <- df_bes%>% 
mutate(partyIdW7.f=fct_recode(partyIdW7,"No-partyId"="No - none" )) 

# We shorten parties names:
  df_bes <- df_bes%>% 
    mutate(partyIdW7.f=fct_recode(partyIdW7.f, "UKIP"="United Kingdom Independence Party (UKIP)"), 
           partyIdW7.f=fct_recode(partyIdW7.f, "SNP"="Scottish National Party (SNP)"), 
           
           partyIdW7.f=fct_recode(partyIdW7.f, "Lib Dem"="Liberal Democrat"), 
           partyIdW7.f=fct_recode(partyIdW7.f, "Conserv."="Conservative"), 
          
           partyIdW7.f=fct_recode(partyIdW7.f, "SNP"="SNP"),
           partyIdW7.f=fct_recode(partyIdW7.f, "Other"="Plaid Cymru"),
           partyIdW7.f=fct_recode(partyIdW7.f, "Other"="Other"),
           partyIdW7.f=fct_recode(partyIdW7.f, "Green P."="Green Party"),
           partyIdW7.f=fct_recode(partyIdW7.f, "Green P."="Green Party"))
    

  table(df_bes$partyIdW7.f)      
  
```


## 1.3. Score varaibles (composite variables)
#### 1.3.1. Media consumption 
```{r}

count_na <- function(x) sum(is.na(x)) 

# This is the function for create the scores. To do this we have calculated the mean. 
score_media<-  df_bes %>% dplyr::select(
infoSourceTVW7,
infoSourcePaperW7,
infoSourceRadioW7,
infoSourceInternetW7) %>% 
  mutate(mean_media = rowMeans(., na.rm = T),
# Those two following lines are used to ignore the "no-answer" of those interviewees who just did not respond tho one item. 
# Ej.: interviewee Nº6 answers: item1=4 , item2=4, item3=NA, item4=1: 
# mean= 3.75
# This interviewee has not responded to item 3 so we have just calculated the mean with the Attitudes of the three other items that were answered. 
# We have proceeded in the same way for the rest of scores varaibles. The underlying reason was not to increase the number of NAs.
         count_media = apply(., 1, count_na)) %>% 
    mutate(score.mediaW7=ifelse(count_media<=1, mean_media, NA))


df_bes$score.mediaW7 <- score_media$score.mediaW7

# To check the distribution of the varaible: 
hist(df_bes$score.mediaW7)

# To check how much NAs has the variable finally. 
summary(df_bes$score.mediaW7)
```




#### 1.3.2. Attitudes score Heath-Evans
```{r}

#######################################
### Authoritarianism ###

# To check the Cronbach's alpha: 
df_bes %>%
dplyr::select(
al1W7_W9,
al2W7_W9,
al3W7_W9,
al4W7_W9,
al5W7_W9) %>% 
  psych::alpha()
# raw_alpha: 0.81

count_na <- function(x) sum(is.na(x)) 

df_al<-  df_bes %>% dplyr::select(
al1W7_W9,
al2W7_W9,
al3W7_W9,
al4W7_W9,
al5W7_W9) %>% 
  mutate(mean_auth = rowMeans(., na.rm = T),
         count_auth = apply(., 1, count_na))%>% 
    mutate(score.auth.W7=ifelse(count_auth<=1, mean_auth, NA))


df_bes$score.auth.W7 <- df_al$score.auth.W7


hist(df_bes$score.auth.W7)

summary(df_bes$score.auth.W7)




############################
### Economic conservatism ###

df_bes %>%
dplyr::select(
lr1W7_W9,
lr2W7_W9,
lr3W7_W9,
lr4W7_W9,
lr5W7_W9) %>% 
  psych::alpha()
# raw_alpha: 0.86

count_na <- function(x) sum(is.na(x)) 

df_lr<-  df_bes %>% 
  mutate(
# With this operation we change the direction of the score to make it fit with the direction of the former value score. 
lr1W7_W9=6-lr1W7_W9,
lr2W7_W9=6-lr2W7_W9,
lr3W7_W9=6-lr3W7_W9,
lr4W7_W9=6-lr4W7_W9,
lr5W7_W9=6-lr5W7_W9
  ) %>% dplyr::select(
lr1W7_W9,
lr2W7_W9,
lr3W7_W9,
lr4W7_W9,
lr5W7_W9) %>% 
  mutate(mean_EcoCon = rowMeans(., na.rm = T),
         count_EcoCon = apply(., 1, count_na))%>% 
    mutate(score.EcoCon.W7=ifelse(count_EcoCon<=1, mean_EcoCon, NA))


df_bes$score.EcoCon.W7 <- df_lr$score.EcoCon.W7


hist(df_bes$score.EcoCon.W7)

summary(df_bes$score.EcoCon.W7)



```


#### 1.3.3. Attitudes score Anti-Egalitarianism
```{r}

#######################################
### Anti-egalitarianism ###

df_bes %>%
dplyr::select(
blackEqualityW6_W14, 
femaleEqualityW6_W14, 
gayEqualityW6_W14
) %>% 
  psych::alpha()
# raw_alpha: 0.79

count_na <- function(x) sum(is.na(x)) 

df_auth_OS<-  df_bes %>% dplyr::select(
blackEqualityW6_W14, 
femaleEqualityW6_W14, 
gayEqualityW6_W14) %>% 
  mutate(mean_AntiEgal = rowMeans(., na.rm = T),
         count_AntiEgal = apply(., 1, count_na))%>% 
    mutate(score.AntiEgal.W7=ifelse(count_AntiEgal<=1, mean_AntiEgal, NA))


df_bes$score.AntiEgal.W7 <- df_auth_OS$score.AntiEgal.W7

hist(df_bes$score.AntiEgal.W7 )
summary(df_bes$score.AntiEgal.W7)

```

#### 1.3.4. Attitudes score Economic conservatism
```{r}

############################
### Economic conservatism ###

df_bes %>%
dplyr::select(
cutsTooFarNationalW7, 
cutsTooFarNHSW7, 
cutsTooFarLocalW7
# privatTooFarW7 # we have deleted this item because in doing so the alpha could be improved. 
) %>% 
  psych::alpha()
# raw_alpha: 0.87 

count_na <- function(x) sum(is.na(x)) 

df_SL_OS<-  df_bes %>% 
  mutate(
cutsTooFarNationalW7=6-cutsTooFarNationalW7,
cutsTooFarNHSW7=6-cutsTooFarNHSW7,
cutsTooFarLocalW7=6-cutsTooFarLocalW7) %>% 
  dplyr::select(
cutsTooFarNationalW7, 
cutsTooFarNHSW7, 
cutsTooFarLocalW7) %>% 
  mutate(mean_EcoCon.Serv = rowMeans(., na.rm = T),
         count_EcoCon.Serv = apply(., 1, count_na))%>% 
    mutate(score.EcoCon.Serv.W7=ifelse(count_EcoCon.Serv<=1, mean_EcoCon.Serv, NA))
 

 
df_bes$score.EcoCon.Serv.W7 <- df_SL_OS$score.EcoCon.Serv.W7

hist(df_bes$score.EcoCon.Serv.W7 )
summary(df_bes$score.EcoCon.Serv.W7)
```



# 2. FINAL DATA FRAME for regressions
```{r}

# We crate a new data frame with the final variables
df_regres<- df_bes %>% select(id,
wt_new_W7,
enviroGrowthW7,
enviroGrowthW17,
enviroProtectionW7,
enviroProtectionW14,

ageW7, 
gender,
educationW7,

p_work_statW7,
econPersonalRetroW7,
econGenRetroW7,
p_gross_householdW7,
occupation, 
subjClassW7.f, 

score.mediaW7,
partyIdW7.f,
leftRightW7,

# Here the principal independent varaibles explained: 
score.auth.W7, # score Attitudes Authoritarian ("auth") 
score.EcoCon.W7, # score Attitudes Economic conservatism ("EcoCon") 

score.AntiEgal.W7, # score Attitudes Authoritarian ("auth") opportunities-services operationalisation ("OS")
score.EcoCon.Serv.W7 # score Attitudes Economic conservatism ("EcoCon.Serv"), services operationalisation ("OS")
)


# We change the order of the Environmental Protection scale dependent varaible. In this way it fits with the direction of the the dependent varaible (Economy  vs Environment). 
# The inversion results in a scale from 1 to 5 in which 1 expresses the opinion "the policies have been excessive" and 5 that "policies have not been sufficient". 

 df_regres<- df_regres %>% 
   mutate(enviroProtectionW7=6-enviroProtectionW7, 
          enviroProtectionW14=6-enviroProtectionW14)
 
 table(df_regres$enviroProtectionW7)
```


## 2.1. Cathegory of reference
```{r}
# We set the cathegory of reference for the factor varaibles. 
 df_regres<- df_regres %>% 
   mutate(gender=relevel(gender, ref = "Male"),
          p_work_statW7=relevel(p_work_statW7, ref = "Unemployed"),
          subjClassW7.f=relevel(subjClassW7.f, ref = "Working-class"),
          partyIdW7.f=relevel(partyIdW7.f, ref = "Labour"), 
          occupation=relevel(occupation, ref = "High-quali"))

```


# 3.  DESCRIPTIVE tables and histograms 

Note: all the tables and plots are referenced to the outputs that can be consulted in the appendix or in the article.

## 3.1. Dependent Variables descriptives
```{r}

# Dependent Variables Table ("T.DV") (see Table A1 in appendix). 
T.DV <- df_regres %>% 
  tableone::CreateTableOne(vars = dplyr::select(df_regres, enviroGrowthW7, enviroGrowthW17,enviroProtectionW7, enviroProtectionW14) %>% names(),
     # strata ="gender", 
    data = .,includeNA = F, # le puedes detallar que te incluya NAs
 test = F, testApprox = chisq.test, # y para incluir un test 
  argsApprox = list(correct = TRUE)) %>%print(varLabels = TRUE)

# For exporting: 
kableExtra:: kbl(T.DV) %>%
  kable_classic(full_width = F, html_font = "Cambria")



#  Dependent Variables Histogram (see Figure A1 in appendix). 
## enviroGrowthW7 Histogram 2016 = w1
h.eg.w7<- df_regres%>% 
ggplot(aes(x=enviroGrowthW7)) +
  geom_histogram(binwidth=0.3,fill="#787878") +
  geom_vline(aes(xintercept=mean(enviroGrowthW7, na.rm = T)),
             color="blue", linetype="dashed", size=0.7)+
labs(title="2016", y="Frequency", x="Economy vs Environment")+
  theme_minimal()
?geom_vline

## enviroGrowthW17 Histogram 2019 = w2
h.eg.w17<- df_regres%>% 
ggplot(aes(x=enviroGrowthW17)) +
  geom_histogram(binwidth=0.3,fill="#5A5A5A") +
  geom_vline(aes(xintercept=mean(enviroGrowthW17, na.rm = T)),
             color="blue", linetype="dashed", size=0.7)+
labs(title="2019", y="", x="Economy vs Environment")+
  theme_minimal()


## enviroProtectionW7 Histogram 2016 = w1
h.ep.w7<- df_regres%>% 
ggplot(aes(x=enviroProtectionW7)) +
  geom_histogram(binwidth=0.3,fill="#787878") +
  geom_vline(aes(xintercept=mean(enviroProtectionW7, na.rm = T)),
             color="blue", linetype="dashed", size=0.7)+
labs(title="2016", y="Frequency", x="Environmental Protection")+
  theme_minimal()

## enviroProtectionW14 Histogram 2018 = w2
h.ep.w14<- df_regres%>% 
ggplot(aes(x=enviroProtectionW14)) +
  geom_histogram(binwidth=0.3,fill="#5A5A5A") +
  geom_vline(aes(xintercept=mean(enviroProtectionW14, na.rm = T)),
             color="blue", linetype="dashed", size=0.7)+
labs(title="2018", y="", x="Environmental Protection")+
  theme_minimal()

# For exporting: 
ggarrange(h.eg.w7, h.eg.w17, h.ep.w7, h.ep.w14, ncol = 2, nrow = 2)

# Export size in PDF
# 5.30 x 5.20

```


## 3.2. Independent Variables descriptives
```{r}
# Independent Variables Table ("T.IV"). See Table A2 in appendix. 

T.VI<-df_regres %>% 
  tableone::CreateTableOne(vars = dplyr::select(df_regres, ageW7,  gender,educationW7,p_work_statW7,subjClassW7.f,p_gross_householdW7,occupation,econPersonalRetroW7,econGenRetroW7,score.mediaW7,partyIdW7.f,leftRightW7,score.EcoCon.W7,score.auth.W7, score.EcoCon.Serv.W7,score.AntiEgal.W7 ) %>% names(),
     # strata ="gender", 
    data = .,includeNA = F, # le puedes detallar que te incluya NAs
 test = F, testApprox = chisq.test, # y para incluir un test 
  argsApprox = list(correct = TRUE)) %>%print(varLabels = TRUE)

# For exporting: 
kableExtra:: kbl(T.VI) %>%
  kable_classic(full_width = F, html_font = "Cambria")
```




## 3.3. Correlations between the different operationalizations of Attitudes
```{r}
# See Table A3 in appendix. 

cor.test(df_regres$score.EcoCon.W7, df_regres$score.EcoCon.Serv.W7)
cor.test(df_regres$score.AntiEgal.W7, df_regres$score.auth.W7)

cor.test(df_regres$score.auth.W7, df_regres$score.EcoCon.Serv.W7)
cor.test(df_regres$score.EcoCon.Serv.W7, df_regres$score.AntiEgal.W7)

```


## 3.4. Intra-person variation
## Cross tables 
The below tables shows interviewers intra individual changes.
#### Dependent Variables
```{r}

Table.EnviroGrowth <- df_bes %>%
  filter(!is.na(enviroGrowthW7)) %>% 
  filter(!is.na(enviroGrowthW17)) %>% 
  tabyl(enviroGrowthW7,enviroGrowthW17) %>% adorn_percentages("row") %>%
  adorn_pct_formatting(digits = 0) 

kableExtra:: kbl(Table.EnviroGrowth) %>%
  kable_classic(full_width = F, html_font = "Cambria")


#########

Table.enviroProtect<- df_regres %>%
  filter(!is.na(enviroProtectionW7)) %>% 
  filter(!is.na(enviroProtectionW14)) %>% 
  tabyl(enviroProtectionW7,enviroProtectionW14) %>% adorn_percentages("row") %>% adorn_pct_formatting(digits = 0) 

kableExtra:: kbl(Table.enviroProtect) %>%
  kable_classic(full_width = F, html_font = "Cambria")


```



# 4. REGRESSION ANALYSIS 

## 4.1. Main regression models
#### 4.1.1. DV: Economy vs Environment  
```{r}



# CROSSECTIONAL MODEL
m.Eco.Env <-lm(enviroGrowthW7 ~ ageW7+  gender+educationW7+p_work_statW7+subjClassW7.f+p_gross_householdW7+occupation+econPersonalRetroW7+econGenRetroW7+score.mediaW7+partyIdW7.f+leftRightW7+score.EcoCon.W7+score.auth.W7+score.AntiEgal.W7                                                 
, data=df_regres, weights = wt_new_W7)


# LAGGED MODEL
m.Eco.Env.lagg <-lm(enviroGrowthW17~ ageW7+gender+educationW7+p_work_statW7+subjClassW7.f+p_gross_householdW7+occupation+econPersonalRetroW7+econGenRetroW7+score.mediaW7+partyIdW7.f +leftRightW7 +score.EcoCon.W7+score.auth.W7+score.AntiEgal.W7                                                 
, data=df_regres, weights = wt_new_W7)


# AUTORREGRESIVE MODEL
m.Eco.Env.autorreg <-lm(enviroGrowthW17~ enviroGrowthW7+ ageW7+gender+educationW7+p_work_statW7+subjClassW7.f+p_gross_householdW7+occupation+econPersonalRetroW7+econGenRetroW7+score.mediaW7+partyIdW7.f +leftRightW7 +score.EcoCon.W7+score.auth.W7+score.AntiEgal.W7                                                 
, data=df_regres, weights = wt_new_W7)


screenreg(list(m.Eco.Env, m.Eco.Env.lagg, m.Eco.Env.autorreg))


# Z-Test to check if the differences between the three attitudes scores are different. 
# Formula: Z=(coef1-coef2)/square root(standar error of coef1^2+standar error of coef2^2)
# Critic values : 1.645 (90%); 1,96 (95%);  2,57 (99%)
# If the result shows numbers above 1,96 this means that the Anti-Egalitarian attitudes significantly have more weight in the models than the other scores. 
# Below are shown in pairs the resoult of each model.  

# CROSSECTIONAL
## AntiEgal VS EcoCon
(-0.33-(-0.42))/sqrt((0.03^2)+(0.03^2))
# Result: 2.12132
## AntiEgal VS auth
(-0.34-(-0.42))/sqrt((0.03^2)+(0.03^2)) 
# Result: 1.885618 
## auth VS EcoCon
(-0.33-(-0.34))/sqrt((0.03^2)+(0.03^2)) 
# Result: 0.2357023 

# LAGGED MODEL
## AntiEgal VS EcoCon
(-0.27-(-0.60))/sqrt((0.05^2)+(0.05^2))
# Result: 4.666905
## AntiEgal VS auth
(-0.25 -(-0.60))/sqrt((0.05^2)+(0.05^2))
# Result: 4.949747
## auth VS EcoCon
(-0.27-(-0.25))/sqrt((0.05^2)+(0.05^2))
# Result: -0.2828427

# AUTORREGRESIVE MODEL
## AntiEgal VS EcoCon
(-0.15-(-0.39))/sqrt((0.04^2)+(0.04^2))
# Result: 4.242641
## AntiEgal VS auth
(-0.12-(-0.39))/sqrt((0.04^2)+(0.04^2))
# Result: 4.772971
## auth VS EcoCon
(-0.15-(-0.12))/sqrt((0.04^2)+(0.04^2))
# Result: -0.5303301



# For exporting the regression table:
  # stargazer(m.Eco.Env, m.Eco.Env.lagg, m.Eco.Env.autorreg, 
          type="html",
          align=TRUE,
          title= "Table A7. Economy  vs Environment,
          style= "ajps",
          out  ="table_A4.html",
          ci.level = 0.9,
          dep.var.labels=c("Economy  vs Environment"),
          covariate.labels=c("Economy vs Environment W7"), 
          single.row=TRUE, 
          digits=2)
 

```



#### 4.1.2. Plots
```{r}


# See Figure 1 in the article.

# Plot with the main results of the section 4.1.1. DV: Economy vs Environment:  
p1<- plot_summs(m.Eco.Env, m.Eco.Env.lagg, m.Eco.Env.autorreg,
           coefs = c("EcoCon." = "score.EcoCon.W7",
                     "Auth." = "score.auth.W7",
                     "AntiEgal." ="score.AntiEgal.W7"), 
           model.names = c("Cross-sectional", "Lagged", "Autoregressive"), x.font.size = 10, 
    colors  = c("red", "black", "grey"),  legend.title = "") + labs(x = "\n Economy vs Environment \n ", y = NULL) 



# Export size in PDF
#  6 x 4,5 
```



#### 4.1.3. DV: Environmental Protection  
```{r}


# CROSSECTIONAL MODEL
m.Env.Prot <-lm(enviroProtectionW7 ~ ageW7+  gender+educationW7+p_work_statW7+subjClassW7.f+p_gross_householdW7+occupation+econPersonalRetroW7+econGenRetroW7+score.mediaW7+partyIdW7.f+leftRightW7+score.EcoCon.W7+score.auth.W7+score.AntiEgal.W7                                                 
, data=df_regres, weights = wt_new_W7)


# LAGGED MODEL
m.Env.Prot.lagg <-lm(enviroProtectionW14 ~ ageW7+gender+educationW7+p_work_statW7+subjClassW7.f+p_gross_householdW7+occupation+econPersonalRetroW7+econGenRetroW7+score.mediaW7+partyIdW7.f +leftRightW7 +score.EcoCon.W7+score.auth.W7+score.AntiEgal.W7                                    
, data=df_regres, weights = wt_new_W7)


# AUTORREGRESIVE MODEL
m.Env.Prot.autorreg <-lm(enviroProtectionW14 ~ enviroProtectionW7+ ageW7+gender+educationW7+p_work_statW7+subjClassW7.f+p_gross_householdW7+occupation+econPersonalRetroW7+econGenRetroW7+score.mediaW7+partyIdW7.f +leftRightW7 +score.EcoCon.W7+score.auth.W7+score.AntiEgal.W7                                                 
, data=df_regres, weights = wt_new_W7)


screenreg(list(m.Env.Prot, m.Env.Prot.lagg, m.Env.Prot.autorreg))


# Z-Test 
# Critic values : 1.645 (90%); 1,96 (95%);  2,57 (99%)

# CROSSECTIONAL
## AntiEgal VS EcoCon
(-0.11-(-0.30))/sqrt((0.01^2)+(0.01^2))
# Result: 13.43503
## AntiEgal VS auth
(-0.08-(-0.30))/sqrt((0.01^2)+(0.01^2))
# Result: 15.55635
## auth VS EcoCon
(-0.11-(-0.08))/sqrt((0.01^2)+(0.01^2))
# Result: -2.12132

# LAGGED MODEL
## AntiEgal VS EcoCon
(-0.11-(-0.27 ))/sqrt((0.02^2)+(0.02^2))
# Result: 5.656854
## AntiEgal VS auth
(-0.11-(-0.27 ))/sqrt((0.02^2)+(0.02^2))
# Result: 5.656854
## auth VS EcoCon
(-0.11-(-0.11 ))/sqrt((0.02^2)+(0.02^2))
# Result: 0

# AUTORREGRESIVE MODEL
## AntiEgal VS EcoCon
(-0.06-(-0.12))/sqrt((0.02^2)+(0.01^2))
# Result: 2.683282
## AntiEgal VS auth
(-0.08-(-0.12))/sqrt((0.02^2)+(0.01^2))
# Result: 1.788854
## auth VS EcoCon
(-0.06 -(-0.08))/sqrt((0.02^2)+(0.02^2))
# Result: 0.7071068

# For exporting the regression table:
  # stargazer(m.Env.Prot, m.Env.Prot.lagg, m.Env.Prot.autorreg, 
          type="html",
          align=TRUE,
          title= "Tabla A8. Environmental Protection",
          style= "ajps",
          out  ="table_A5.html",
          ci.level = 0.9,
          dep.var.labels=c("Environmental Protection"),
          covariate.labels=c("Environmental Protection W7"), 
          single.row=TRUE, 
          digits=2)


```



#### 4.1.4. Plots
```{r}

# See Figure 1 in the article.



# Plot with the main results of the section 4.1.2. DV: Environmental Protection:  
p2<-  plot_summs(m.Env.Prot, m.Env.Prot.lagg, m.Env.Prot.autorreg,
           coefs = c("EcoCon." = "score.EcoCon.W7",
                     "Auth." = "score.auth.W7",
                     "AntiEgal." = "score.AntiEgal.W7"), 
           model.names = c("Cross-sectional", "Lagged", "Autoregressive"), x.font.size = 15, 
    colors  = c("red", "black", "grey"),  legend.title = "") + labs(x = "\n Environmental Protection \n ", y = NULL) 



ggarrange(p1 + font("title", size = 13)+font("x.text", size = 10),
          
          p2 + font("title", size = 13)+font("x.text", size = 10),
          common.legend=T, legend=c("bottom"))  


Export size in PDF
#  6 x 4,5 
```




## 4.2. Robustness check 
Here we run the same models but using a different operationalisation of Economic Conservatism. 

#### 4.2.1. DV: Economy vs Environment  
```{r}

# CROSSECTIONAL MODEL
m.Eco.Env.serv <-lm(enviroGrowthW7 ~ ageW7+  gender+educationW7+p_work_statW7+subjClassW7.f+p_gross_householdW7+occupation+econPersonalRetroW7+econGenRetroW7+score.mediaW7+partyIdW7.f+leftRightW7+score.EcoCon.Serv.W7+score.auth.W7+score.AntiEgal.W7                                                 
, data=df_regres, weights = wt_new_W7)


# LAGGED MODEL
m.Eco.Env.lagg.serv  <-lm(enviroGrowthW17~ ageW7+gender+educationW7+p_work_statW7+subjClassW7.f+p_gross_householdW7+occupation+econPersonalRetroW7+econGenRetroW7+score.mediaW7+partyIdW7.f +leftRightW7 +score.EcoCon.Serv.W7+score.auth.W7+score.AntiEgal.W7                                                 
, data=df_regres, weights = wt_new_W7)


# AUTORREGRESIVE MODEL
m.Eco.Env.autorreg.serv  <-lm(enviroGrowthW17~ enviroGrowthW7+ ageW7+gender+educationW7+p_work_statW7+subjClassW7.f+p_gross_householdW7+occupation+econPersonalRetroW7+econGenRetroW7+score.mediaW7+partyIdW7.f +leftRightW7 +score.EcoCon.Serv.W7+score.auth.W7+score.AntiEgal.W7                                                 
, data=df_regres, weights = wt_new_W7)


screenreg(list(m.Eco.Env.serv, m.Eco.Env.lagg.serv, m.Eco.Env.autorreg.serv))


# Z-Test to check if the differences between the three attitudes scores are different. 
# Formula: Z=(coef1-coef2)/square root(standar error of coef1^2+standar error of coef2^2)
# Critic values : 1.645 (90%); 1,96 (95%);  2,57 (99%)
# If the result shows numbers above 1,96 this means that the Anti-Egalitarian attitudes significantly have more weight in the models than the other scores. 
# Below are shown in pairs the resoult of each model.  

# CROSSECTIONAL
## AntiEgal VS EcoCon
(-0.29-(-0.42))/sqrt((0.03^2)+(0.03^2))
# Result: 3.064129
## AntiEgal VS auth
(-0.30-(-0.43))/sqrt((0.03^2)+(0.03^2)) 
# Result: 3.064129 
## auth VS EcoCon
(-0.29 -(-0.30))/sqrt((0.03^2)+(0.03^2)) 
# Result: 0.2357023 

# LAGGED MODEL
## AntiEgal VS EcoCon
(-0.34-(-0.59))/sqrt((0.05^2)+(0.05^2))
# Result: 3.535534
## AntiEgal VS auth
(-0.25 -(-0.59))/sqrt((0.05^2)+(0.05^2))
# Result: 4.808326
## auth VS EcoCon
(-0.34 -(-0.25))/sqrt((0.05^2)+(0.05^2))
# Result: -1.272792

# AUTORREGRESIVE MODEL
(-0.23-(-0.38))/sqrt((0.04^2)+(0.04^2))
# Result: 2.65165
(-0.12-(-0.38))/sqrt((0.04^2)+(0.04^2))
# Result: 4.596194
## auth VS EcoCon
(-0.23-(-0.12))/sqrt((0.04^2)+(0.04^2))
# Result: -1.944544



# For exporting the regression table:
    # stargazer(m.Eco.Env.serv, m.Eco.Env.lagg.serv, m.Eco.Env.autorreg.serv, 
          type="html",
          align=TRUE,
          title= "Table A9. Economy  vs Environment – Heath-Evans",
          style= "ajps",
          out  ="table_A6.html",
          ci.level = 0.9,
          dep.var.labels=c("Economy  vs Environment"),
          covariate.labels=c("Economy vs Environment W7"), 
          single.row=TRUE, 
          digits=2)
 

```



#### 4.2.2. Plots
```{r}


# See Figure 1 in the article.

# Plot with the main results of the section 4.1.1. DV: Economy vs Environment:  
p3 <- plot_summs(m.Eco.Env.serv, m.Eco.Env.lagg.serv, m.Eco.Env.autorreg.serv,
           coefs = c("EcoCon." = "score.EcoCon.Serv.W7",
                     "Auth." = "score.auth.W7",
                     "AntiEgal." ="score.AntiEgal.W7"), 
           model.names = c("Cross-sectional", "Lagged", "Autoregressive"), x.font.size = 10, 
    colors  = c("red", "black", "grey"),  legend.title = "") + labs(x = "\n Economy vs Environment \n ", y = NULL) 



# Export size in PDF
#  6 x 4,5 
```



#### 4.2.3. DV: Environmental Protection  
```{r}



# CROSSECTIONAL MODEL
m.Env.Prot.serv <-lm(enviroProtectionW7 ~ ageW7+  gender+educationW7+p_work_statW7+subjClassW7.f+p_gross_householdW7+occupation+econPersonalRetroW7+econGenRetroW7+score.mediaW7+partyIdW7.f+leftRightW7+score.EcoCon.Serv.W7+score.auth.W7+score.AntiEgal.W7                                                 
, data=df_regres, weights = wt_new_W7)


# LAGGED MODEL
m.Env.Prot.lagg.serv <-lm(enviroProtectionW14 ~ ageW7+gender+educationW7+p_work_statW7+subjClassW7.f+p_gross_householdW7+occupation+econPersonalRetroW7+econGenRetroW7+score.mediaW7+partyIdW7.f +leftRightW7 +score.EcoCon.Serv.W7+score.auth.W7+score.AntiEgal.W7                                    
, data=df_regres, weights = wt_new_W7)


# AUTORREGRESIVE MODEL
m.Env.Prot.autorreg.serv <-lm(enviroProtectionW14 ~ enviroProtectionW7+ ageW7+gender+educationW7+p_work_statW7+subjClassW7.f+p_gross_householdW7+occupation+econPersonalRetroW7+econGenRetroW7+score.mediaW7+partyIdW7.f +leftRightW7 +score.EcoCon.Serv.W7+score.auth.W7+score.AntiEgal.W7                                                 
, data=df_regres, weights = wt_new_W7)


screenreg(list(m.Env.Prot.serv, m.Env.Prot.lagg.serv, m.Env.Prot.autorreg.serv))


# Z-Test 
# Critic values : 1.645 (90%); 1,96 (95%);  2,57 (99%)

# CROSSECTIONAL
## AntiEgal VS EcoCon
(-0.10-(-0.30))/sqrt((0.01^2)+(0.01^2))
# Result: 14.14214
## AntiEgal VS auth
(-0.07-(-0.30))/sqrt((0.01^2)+(0.01^2))
# Result: 16.26346
## auth VS EcoCon
(-0.10-(-0.07))/sqrt((0.01^2)+(0.01^2))
# Result: -2.12132

# LAGGED MODEL
## AntiEgal VS EcoCon
(-0.05-(-0.27 ))/sqrt((0.02^2)+(0.02^2))
# Result: 7.778175
## AntiEgal VS auth
(-0.10-(-0.27 ))/sqrt((0.02^2)+(0.02^2))
# Result: 6.010408
## auth VS EcoCon
(-0.05 -(-0.10 ))/sqrt((0.02^2)+(0.02^2))
# Result: 1.767767

# AUTORREGRESIVE MODEL
## AntiEgal VS EcoCon
(-0.00-(-0.12))/sqrt((0.01^2)+(0.01^2))
# Result: 8.485281
## AntiEgal VS auth
(-0.08-(-0.12))/sqrt((0.02^2)+(0.01^2))
# Result: 1.788854
## auth VS EcoCon
(-0.00-(-0.08))/sqrt((0.01^2)+(0.02^2))
# Result: 3.577709

# For exporting the regression table:
   # stargazer(m.Env.Prot.serv, m.Env.Prot.lagg.serv, m.Env.Prot.autorreg.serv, 
          type="html",
          align=TRUE,
          title= "Tabla A10. Environmental Protection - Heath-Evans",
          style= "ajps",
          out  ="table_A7.html",
          ci.level = 0.9,
          dep.var.labels=c("Environmental Protection"),
          covariate.labels=c("Environmental Protection W7"), 
          single.row=TRUE, 
          digits=2)


```



#### 4.2.4. Plots
```{r}

# See Figure 1 in the article.



# Plot with the main results of the section 4.1.2. DV: Environmental Protection:  
p4 <- plot_summs(m.Env.Prot.serv, m.Env.Prot.lagg.serv, m.Env.Prot.autorreg.serv,
           coefs = c("EcoCon." = "score.EcoCon.Serv.W7",
                     "Auth." = "score.auth.W7",
                     "AntiEgal." = "score.AntiEgal.W7"), 
           model.names = c("Cross-sectional", "Lagged", "Autoregressive"), x.font.size = 15, 
    colors  = c("red", "black", "grey"),  legend.title = "") + labs(x = "\n Environmental Protection \n ", y = NULL) 


ggarrange(p3 + font("title", size = 13)+font("x.text", size = 10),
          
          p4 + font("title", size = 13)+font("x.text", size = 10),
          common.legend=T, legend=c("bottom"))  



# Export size in PDF
#  6 x 4,5 
```

